The Trapped Polarized Fermi Gas at Unitarity 
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We consider population-imbalanced two-component Fermi gases under external harmonic confine- 
ment interacting through short-range two-body potentials with diverging s-wave scattering length. 
Using the fixed-node diffusion Monte Carlo method, the energies of the "normal state" are deter- 
mined as functions of the population-imbalance and the number of particles. The energies of the 
trapped system follow, to a good approximation, a universal curve even for fairly small systems. A 
simple parameterization of the universal curve is presented and related to the equation of state of 
the bulk system. 
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Strongly-correlated few- and many-body systems play 
an important role in atomic, nuclear and condensed mat- 
ter physics (see, e.g., Refs. [EH; IS El)- These systems are 
characterized by intricate particle-particle correlations, 
which can be difficult to capture by mean-field frame- 
works. Furthermore, the development of effective beyond 
mean-field approaches is complicated by the fact that a 
small parameter can typically not be identified. Thus, 
the most promising avenues for the theoretical descrip- 
tion of strongly-correlated systems start with a micro- 
scopic many-body Hamiltonian. In a few fortuitous cases, 
such as strongly-correlated one-dimensional systems [Bj], 
exact analytical solutions can be found for certain classes 
of model Hamiltonian. In other cases, however, a quan- 
titative description of strongly-correlated systems relies 
on numerical approaches. 

This paper treats strongly-correlated two-component 
Fermi gases under external harmonic confinement with 
diverging s-wave scattering length a s at zero temperature 
using a numerical Monte Carlo approach. The masses of 
the two species are assumed to be the same, and the 
properties of the system are determined as a function of 
the population difference between the two components. 
The interest in population-imbalanced atomic Fermi sys- 
tems [I, 0, H, @, USi El 03 is inherently linked to the 
mismatch of the Fermi surfaces of the two components, 
which has, e.g., been predicted to lead under certain cir- 
cumstances to pairing with non-zero momentum, i.e., to 
the formation of so-called FFLO states [HI, [3] • So far, 
the experimental search for FFLO states has, however, 
not been met with success. 

Population-imbalanced Fermi systems have been real- 
ized using ultracold atom samples, and a rich behavior as 
functions of population-imbalance and temperature has 
been observed [Tl|, [l2j ■ In these experiments, composite 
fermionic atoms such as 6 Li are trapped in two different 
internal hyperfine states, referred to as spin-up and spin- 
down atoms in the following. The population-imbalance 
or polarization can be adjusted straightforwardly, and 
the interspecies scattering length a s between the spin- 
up and spin-down fcrmions can be varied by applying 
an external magnetic field in the vicinity of a so-called 
Fano-Feshbach resonance. To date, all experiments on 



population-imbalanced Fermi gases have been performed 
with comparatively large atom samples. Here, we con- 
sider the properties of small systems, which can be real- 
ized using present-day technology by loading a degener- 
ate Fermi gas into a deep optical lattice, for which the 
tunneling between neighboring sites is neligible. 

We determine the energetics of small trapped two- 
component Fermi gases at unitarity with varying 
population-imbalance using the fixed-node diffusion 
Monte Carlo (FN-DMC) method [Hj]. Two different pa- 



rameterizations of the many-body nodal surface are con- 
sidered: The nodal surface of the non-interacting trapped 
gas, and a nodal surface that is constructed by anti- 
symmetrizing a set of pair functions. Not surprisingly, 
the trapped gas described by the former nodal surface 
can be related to the normal state of the homogeneous 
system. We find that the energies of the "trapped normal 
gas" with varying population-imbalance and number of 
particles fall on a universal curve that can be parame- 
terized quite well by three parameters. Remarkedly, the 
energies for systems with as few as N = 5 fermions fall 
on the same universal curve as those for larger systems 
(the largest N considered in this work is 55). The rela- 
tionship between the trapped and homogeneous systems 
is also analyzed by considering structural properties. For 
small N, we find that the energies of the trapped normal 
system are lower than those obtained for a nodal surface 
that accounts for pairing physics. For larger N, in con- 
trast, the nodal surface that accounts for pairing physics 
results in a lower energy for small \N± — iV^I, where Ni 
and N2 denote the number of spin-up and spin-down 
fcrmions, respectively. Our ab initio results for trapped 
population-imbalanced Fermi gases may serve as bench- 
marks for other numerically less involved techniques and 
aid in testing phenomenological models. 

Our model Hamiltonian for the trapped two- 
component Fermi gas with N mass m fermions, where 
N = ^+N 2 and Ni > N 2 , reads 
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Here, r, denotes the position vector of the ith atom and 
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u> the angular trapping frequency. The interaction poten- 
tial V t b between spin-up and spin-down atoms depends on 
the interparticle distance m, =\ri — Tj\, and is char- 
acterized by the interspecies s-wave scattering length a s . 
We model Vtb by a square well potential with range i?o 
and depth V (Vo > 0), Vtb(r) = —Vo for r < R and 
for r > Rq. To describe the unitary regime, we fix Ro 
at Rq = O.Ola/jo, where a^o denotes the harmonic oscil- 
lator length, aho = y/h/ (mu>), and adjust the depth Vq 
so that the free-space two-particle system supports a sin- 
gle zero-energy s-wave bound state. We consider regimes 
away from an intra-species p-wave Feshbach resonance, 
and treat like atoms as non-interacting. 

The ground state wave function for fermions is, as 
a consequence of the Pauli exclusion principle, charac- 
terized by a complicated nodal surface, i.e., the many- 
body wave function changes sign when either two spin-up 
atoms or two spin-down atoms are exchanged. This sign 
change often times leads to inefficient sampling schemes, 
a phenomenon commonly referred to as "fermionic sign 
problem" . To avoid the sign problem, we adopt the FN- 
DMC method [lq . which determines the eigenenergy of 
a state that has the same nodal structure as tpx but 
that may differ from ?/>t in other regions of configuration 
space. In most applications of the FN-DMC method, 
great effort is placed on optimizing ipx so as to obtain a 
tighter upper bound to the true eigenenergy. Our calcu- 
lations that utilize a pairing function fall into this cate- 
gory. Our study of the trapped normal state, in contrast, 
considers the properties of the many-body system for a 
fixed, non-optimized nodal structure. 

The guiding function ipTi of the trapped normal gas 
is written in terms of the ground state wave function 
ipffi of the non-interacting trapped Fermi gas and three 
Jastrow factors Jkk'i ipri — ipNiJnJ22Ji2- The posi- 
tive definite Jastrow factors Jkk' account for correlations 
between atoms from component k and atoms from com- 
ponent k'; they reduce the statistical uncertainties but 
do not alter the nodal surface of ipTi- The function ipNi 
can be written as a product of two Slater determinants, 
ipNi = Det(Mi) x Det(Af2), where the r/th element of 
the matrix Mj,, k = 1 or 2, is given by the single particle 
harmonic oscillator function 4>i(rj) with i,j = 1, • • • , Nk- 
The subscript i labels the excitations: i = 1 corresponds 
to the ground state, i = 2 to the first excited state, and 
so on. For closed shells, i.e., for Nk = 1, 4, 10, 20, 35, • • • , 
the Slater determinant Det(Mfc) is uniquely defined. For 
open shells, in contrast, degenerate states exist and tpNi 
is not uniquely defined. In this case, we consider deter- 
minants build from different sets of harmonic oscillator 
orbitals and report the lowest energy. For Nk = 5, e.g., 
we fill the first two shells and place the remaining parti- 
cle in either a / = or I = 2 orbital (/ denotes the orbital 
angular momentum of <f>i). 

Tables U and [IT] summarize our FN-DMC energies 
for the guiding function tpTi- In addition to the 
FN-DMC energy, we calculate the expectation value 
Evmc, E V mc = {ipT\H\ip T )/(ip T \ipT), of the many- 



TABLE I: FN-DMC energies in units of hcj for two- 
component unitary Fermi gases with various Ni , N2 combi- 
nations (N < 20 and | JVi — Ns\ < 10) calculated using V'ti- 
The energies are uncertain in the last digit reported. The en- 
ergy of the (3, 1) system is 6.60hu; the energies of the other 
small systems (JVi, JVb < 3) are reported in Refs. [l6l. \l7§. 
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31.2 


33.1 
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41.2 


43.2 



TABLE II: FN-DMC energies in units of hu for unitary Fermi 
gases with closed shells and N > 20, calculated using ipTi- 
The energies are uncertain in the last digit reported. 
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body Hamiltonian using the variational Monte Carlo 
(VMC) method [l5j]. For all systems considered in Ta- 
ble HJ [as well as for the (N U N 2 ) = (20,1), (35,1), (20,4), 
and (35,4) systems], EvmC is, for the simulation times 
employed, positive and less than about 15 % higher than 
the corresponding FN-DMC energy. For the systems with 
(N 1: N 2 ) = (20,10), (20,20), (35,10) and (35, 20), in con- 
trast, the VMC energy becomes negative after a fairly 
small number of sampling steps despite the fact that 
the initial configurations correspond to a gas-like sys- 
tem. To obtain the FN-DMC energies for these systems, 
we start the FN-DMC calculations from non-equilibrated 
gas- like configurations and not, as done for the smaller 
systems, from configurations that are distributed accord- 
ing to ipTV T ne resulting FN-DMC energies correspond 
to a gas-like state and appear converged. For the (20, 10) 
system, e.g., we have checked that different initial config- 
urations result, within errorbars, in the same FN-DMC 
energy. The existence of many-body bound states with 
negative energy for large systems interacting through 
finite-ran ge t wo-body potentials is not surprising (see, 
e.g., Ref. Pi|)- In fact, even some of the smaller systems 
with N > N cr may possess a tightly-bound molecular- 
like ground state that is not sampled in our simulations. 
For systems with Ni = N2, it has been found previously 
that the critical number N cr is larger than 6 |lq|: for 
systems with Ni — N2 = 1, it is larger than 5 [lq . 

Figure EQ shows the FN-DMC energies E from Tables UJ 
andlTTl scaled by the total energy Eni of the correspond- 
ing non-interacting gas, as a function of the concentra- 
tion x, where x — N2/N1. A fully polarized, single- 
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FIG. 1: (Color online) The scaled energies E/Eni for the 
trapped normal Fermi gas at unitarity, shown as a function 
of the concentration x, fall to a good approximation on a uni- 
versal curve. The inset shows E/Eni for x = 1 as a function 
of 1/iVi, JVi = 5 - 20; the LDA prediction = 0.75 is 

shown by a solid line. 



FIG. 2: (Color online) Energetics of the trapped population- 
imbalanced normal Fermi gas. Main figure: Symbols show the 
quantity (E — Eni) / Eni,i as a function of x. Inset: Symbols 
show (E — Eni) I 'Eni ,i for a single down-fermion immersed 
in a cloud of up-fermions as a function of 1/iVi (Ni = 5 — 35). 
The dashed lines show Eq. (J2j) with A» = 0.72. The solid line 
shows Eq. © with A* = 0.72, C* = 0.16 and £>* = 1.7. 



component Fermi gas corresponds to x = while a Fermi 
gas with equal number of particles in the two compo- 
nents corresponds to x = 1. Notably, the scaled ener- 
gies E/Eni fall to a good approximation on a universal 
curve. Applying the local density approximation (LDA) 
to the equation of state of the population-balanced ho- 
mogeneous normal system at unitarity, the energy of 
the trapped unitary system with x = 1 is given by 
y/^Epfi, where £ = 0.54 [l9| [a more recent study re- 
ports £ = 0.56(1) Q, implying V£ = 0.75(2)]. The inset 
of Fig. [1] shows the scaled energies E/Eni for x = 1 at 
unitarity as a function of 1/N\ (N\ — 5 — 20). The fact 
that E/Eni is close to 0.75 suggests that the trapped 
population-balanced unitary gas, described by a wave 
function whose nodal surface coincides with that of the 
non-interacting gas, can be described quite accurately 
by applying the LDA to the normal state of the homoge- 
neous system. 

To analyze the energies of trapped polarized unitary 
Fermi gases further, we consider the energy of a single 
down-fermion or impurity immersed in a cloud of up- 
fermions. Symbols in the inset of Fig. [2] show the scaled 
energy difference (E — Eni)/Eni,i, where E is the FN- 
DMC energy of the (Ni, 1) system (Ni = 5 - 35), as a 
function of 1/Ni; since N2 — 1, the quantity l/Ni coin- 
cides with the concentration x. ENi,k denotes the energy 
of the non-interacting single-component Fermi gas with 
Nk atoms (Eni = En 1,1 + Eni^)- If the system was 
non-interacting, the quantity E — Eni would be zero. 
The deviation from zero can thus be interpreted as the 
attractive interaction or binding energy between the im- 
purity and the up-fermions. To determine this energy, 
we write 

(E-E NI )/Em,i = -A*/N 1 . (2) 



A linear fit to the FN-DMC energies, shown by a dashed 
line in the inset of Fig. [2l gives A* = 0.72(1). The anal- 
ysis performed here for the trapped system is similar to 
that performed for the homogeneous system 0, H, H|. 
However, while we fix the impurity mass at its "bare" 
value, work on the homogeneous system treats the im- 
purity as a quasi-particle with effective mass m e ff and 
finds m e ff — 1.09(2)m at unitarity 0]. 

Symbols in the main part of Fig. [5] show the scaled 
energy (E — Eni) /En 1,1 as a function of x (the plot in- 
cludes all energies reported in Tables HI and ILT]) . A dashed 
line shows the quantity — A^x. For larger x, the scaled 
energies lie above the dashed line, indicating a shield- 
ing of the attractive up-down interaction, or equivalently, 
the presence of an effective repulsion between the down- 
fermions. A simple expression for the energy of the nor- 
mal state of trapped two-component systems with arbi- 
trary x (x > l/Aq) reads 

(E - E N i)/Eni,i = -A*x + C*x D '. (3) 

Using A* = 0.72, a fit to our data gives C* = 0.16(1) 
and = 1.7(1). This fit (solid line in Fig. [5]) provides a 
good description of our numerically determined energies. 
Note that our D* value would be somewhat larger if we 
enforced Eq. ([3]) to reproduce the LDA value of 0.75 at 
x= 1. 

In addition to the energies, we analyze the structural 
properties of population-imbalanced Fermi gases. The 
LDA predicts to leading order in x that the density of the 
up-fermions is unchanged while the down-fermions feel 
a modified trapping potential with effective angular fre- 
quency uj eff = (1 + 3A hom /5)m/m e ff 0, [l0| , where 
Ahom = 0.99(1) 9]. The single down-fermion immersed 
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in a cloud of up-fermions is thus predicted to be de- 
scribed by the ground state harmonic oscillator function 
with width aho, e ff, where a ho ,eff = \A/( m e//We//)- A 
fit of the square of the harmonic oscillator function to 
the density profiles of the down-fermion for the (10, 1) 
and (20, 1) systems, calculated using the mixed estima- 
tor [15| , provides a good description of the density profiles 
and gives ciho,eff = 0.88a; lo and 0.86ah o ,ef f, respectively, 
compared with the LDA prediction of a/, ,e// = 0.87ah o . 
The good agreement suggests that the LDA provides a 
valid description even of fairly small trapped systems. 

In addition to the guiding function ipTi, we con- 
sider the guiding ipT2, which is constructed by anti- 
symmetrizing a set of N 2 pair functions / and Ni — N 2 
mutually orthogonal harmonic oscillator functions (re- 
call, Ni > n 2 ) pi mj. It has been shown previ- 
ously [H, [ItJ that ipT2 results in a lower FN-DMC energy 
than tpTi for population-balanced systems with N > 6 
and for odd N systems with \Ni - N 2 \ = 1 and N > 13. 
Here, we extend the analysis to systems with larger 
N% — N 2 . For the (7,5) system, tpT2 and V"ri result, 
within errorbars, in the same energy. For the next larger 
systems with Ni — N 2 = 2, the guiding function ip T2 
results in lower FN-DMC energies than tpxi [E = 28.6, 
33.7, and 38.9haj for the (8, 6), (9, 7) and (10, 8) systems, 



respectively]. For N 1 — N 2 > 3 and N < 20, in contrast, 
we find that the guiding function ip^i results in lower 
FN-DMC energies than vb T2 . 

In summary, this paper treats trapped polarized 
two-component Fermi gases interacting through short- 
range two-body potentials at unitarity. This strongly- 
correlated regime has attracted a great deal of attention 
since the only meaningful length scale is the system's size 
(see, e.g., Ref. 20]); consequently, the properties of the 
gas are governed by a few universal parameters. Using 
the FN-DMC method, we have determined the energies 
of the trapped normal gas as functions of N and x. Our 
energies fall on a universal curve that is well described 
by three parameters. Guiding functions that account for 
pairing physics are also considered and found to result in 
lower energies than those obtained for the trapped nor- 
mal state when N is sufficiently large and Ni — N 2 suffi- 
ciently small. Our results may aid in assessing the accu- 
racy of other numerical approaches such as density func- 
tional theory approaches 21|. Furthermore, they may 
guide optical lattice experiments on ultracold fermionic 
gases. 
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